A frozen glass phase in the multi-index matching problem 
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The multi-index matching is an NP-hard combinatorial optimization problem; for two indices it 
reduces to the well understood bipartite matching problem that belongs to the polynomial com- 
plexity class. We use the cavity method to solve the thermodynamics of the multi-index system 
with random costs. The phase diagram is much richer than for the case of the bipartite matching 
problem: it shows a finite temperature phase transition to a completely frozen glass phase, similar 
to what happens in the random energy model. We derive the critical temperature, the ground state 
energy density, and properties of the energy landscape, and compare the results to numerical studies 
based on exact analysis of small systems. 

PACS numbers: 75.10.Nr, 75.40.-s, 75.40.Mg 



It has been recognized early on that one important mo- 
tivation of the research in spin glass theory is the ubiquity 
of systems with frustration and disorder JJ. In particu- 
lar, recent statistical physics studies have brought inter- 
esting new results in some important computer science 
problems. Notable examples are found in optimization 
(e.g. K-satisfiabihty (/iT-SAT) graph coloring 0, 
or vertex cover M) and information theory (e.g. error 
correcting codes |5j). So far, the most interesting appli- 
cations of spin glass theory have been obtained in this 
emerging field, which witnesses an upsurge of interdisci- 
plinary studies involving physicists, computer scientists, 
and probabilists. 

One of the first optimization problems studied ana- 
lytically by physics methods was the random Bipartite 
Matching Problem (BMP). It is also a simple problem: 
from the computer science point of view, it belongs to 
the class P of polynomial complexity; from the physics 
point of view, it has no phase transition at finite temper- 
ature, and its solution with the replica method shows 
a simple replica symmetric behavior. Interestingly, the 
validity of this solution has been recently confirmed by a 
rigorous mathematical study 0. 

In this work we study the Multi-Index Matching Prob- 
lem (MIMP), a natural extension [1] of the BMP. This is 
a more complicated problem: it belongs to the NP-hard 
class, and as we will see it also exhibits a finite tem- 
perature phase transition, with a low temperature glassy 
phase. Using the cavity method, we find that this phase 
consists of isolated configurations, and we conjecture that 
our method yields an exact solution to this problem. Be- 
cause of its structural resemblance to the BMP, one may 
hope that the MIMP will also be amenable to rigorous 
study, generalizing the construction of to a problem 
with a glass phase. 

The random MIMP — In the BMP one is given two 
sets of M points, and ^2. Each point of 5*1 must be 
"matched" or assigned to one point of 5*2. This matching 
must be one-to-one, and it can be represented by the 



"occupation" of the edges between the points of the two 
sets; we define JT-ij.ij = 1 if the points (ii, 12) G SiX S2 are 
matched, while fM-^^i^ — otherwise. To each matching 
we associate a cost or energy, which is the sum of the 
costs of each occupied edge. 

The MIMP is a straightforward generalization of this 
problem to more than two indices. Given d sets Si,. . . ,Sd, 
each of M sites, a hyperedge is a d-uplet where exactly 
one site from each set appears. For each hyperedge we 



introduce a cost ij 



and the total cost of a (multi- 



index) matching is given, in terms of the occupation num- 
bers of hyperedges, by: 



,id"'ti,...,id ■ 



(1) 



The occupation numbers of hyperedges, ni^,,,,^i^ G {0, 1} 
must be such that each site appears exactly once: 



Vre[l,(i], 



E 



= 1. 



(2) 



The optimization problem consists in finding the min- 
imum cost matching. What makes this problem difficult 
is the constraint jSJ of having each site appear exactly in 
just one hyperedge; for d > 3 the MIMP is NP-hard 
MIMP arise for instance when assigning tasks (jobs) to 
people in particular time slots or in different places. An 
application also arises in the context of track reconstruc- 
tion |loj : given the positions of M unlabeled particles at 
d different times, one is to determine the tracks or tra- 
jectories of each. This kind of formulation is in fact used 
in track reconstruction in high energy physics . 

We shall study the random MIMP where the individ- 
ual costs ii-^.,,,id are independent identically distributed 
random variables. For definiteness we shall take ^^^^...i^ 
to have uniform distribution in [0, 1], although other dis- 
tributions can be studied similarly. 

For a given sample characterized by the values of 
the partition function at inverse temperature /3 = 



2 



is 



-/3ff[{n.,,...,.J] 



(3) 



where /3 is the inverse temperature and the sum is over aU 
possible matchings. In the thermodynamic hmit (M 
oo), only the behavior at the lowest values of £i-^^,,,^i^ 
matters. Indeed, if we consider a given site in any of the d 
sets, it is to be assigned to a low cost hyperedge; generally 
it is possible to assign it to one of the first shortest such 
hyperedges. This means that at large M, the typical 
cost of an occupied hyperedge in the low temperature 
regime should scale as 1/M'^~^. It is thus convenient 
to work with rescaled quantities that are extensive (i.e. 
proportional to M): 



E = W^-^H . 



(4) 



This amounts to considering thermodynamic quantities 
and having (3 = P/M'^~^ as the control parameter: one 
should keep P fixed when taking the large M limit. 

Given these considerations, we conjecture that the free 
energy density is self-averaging as in most disordered sys- 
tems, and in particular as rigorously proved for d = 2 

Numerical study of the ground state — For a given 
sample of the quenched disorder, we want to determine 
the ground state energy Eq which is the minimum of all 
An exhaustive search over all matchings 
works only for very small M (typically Af < 6 when 
d > 3) because of the rapid growth of the number of 
legal matchings, in (M!)''"'^. We have followed instead 
a branch and bound (B&B) approach which allows us to 
study intermediate M. From such an algorithm, we can 
test numerically whether Eq is self-averaging and study 
its large M limit. 

The determination of the best matching uses a search 
tree. All the nodes at level p of this tree correspond to 
having chosen hyperedges for the first p points of the 
set Si (ordered arbitrarily) . To go from level p to level 
{p + 1), we branch on all possible M'^~^ choices for the 
next hyperedge. Then a path from the tree's root (level 
0) to a leaf (level M) is a choice of M hyperedges which 
may or not correspond to a legal matching. The cost of 
a node in the tree is defined as: the sum of the costs of 
its associated hyperedges when they don't overlap, or oo 
if the hyperedges overlap (i.e. they use a point of the d 
sets Si more than once) . 

The B&B algorithm searches the tree and implements 
pruning. For this, it needs an upper bound E^h on Eq; 
we initialize this quantity before performing the search 
via the cost of a legal matching obtained by a greedy as- 
signment of the hyperedges. Then the algorithm starts 
at the root of the tree (level 0) and searches it recur- 
sively. At each level, one branches on the M'^~^ choices 
of hyperedges that take one to the next level. Every time 
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FIG. 1: Mean ground state energy density Eo/M as a func- 
tion of 1/M in the 3- index problem and our best fit. Inset: 
standard deviation (t(_Bo)/v'a7 versus 1/M. 



the current node has a cost greater than i?ub, all of its 
descendent nodes can be ignored as they cannot contain 
the ground state. If we reach level M, we have a legal 
matching which we keep if its cost is less than E^]^ (and 
we update _Eub accordingly) . Upon termination, we have 
the ground state and Eq = i?ub- 

We have implemented this algorithm along with a num- 
ber of optimizations. For our computer program, one 
d — 3 sample at M = 20 takes typically 5 seconds on a 
2 GHz PC, and the CPU time grows by a factor around 
2.2 every time M is increased by 1. We have performed 
runs for M < 22 with 20000 samples at each M. From 
these data, we have extracted Eq , the disorder average of 
Eq; the mean cost per node is shown in Fig.^ The data 
for Af > 10 are well fitted by a quadratic curve in 1/M, 
giving Eq/M — > 3.06 ± 0.03; a power law fit of the same 
quality gives 'E^/M 3.09 ± 0.03. 

In the inset of the figure, we show that the standard 
deviation a{EQ/M) scales as 1/^/M, which is evidence for 
self-averaging and also suggests a central limit theorem 
behavior when M ^ oo. 

Finally, we have also investigated a bit the case of d = 
4; however, we were limited to Af < 15 and used only 
5000 samples. (The CPU time grows by about the same 
factor when Af is increased by 1 as when d = 3.) Our 
best fit in this case leads to Eq/M 7.2(3). 

Thermodynamics and the cavity approach — The re- 
cent formulation of the cavity method 0| for diluted sys- 
tems offers a choice tool to study the thermodynamics of 
the MIMP analytically. Building on the idea that the op- 
timal matching selects preferentially the hyperedges with 
the lowest costs, we dilute the initially complete hyper- 
graph by suppressing hyperedges with £ii^...Sd > CAf ^"'^ 
|l4|. In the resulting graph, the degree of each site is a 
Poisson distributed variable of mean C. When increas- 
ing Af to Af + 1, a new serie of d sites is added. Each 
of them is connected to a finite number of neighbours. 
The partition function of one new site is easily computed 
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FIG. 2: Free energy density as a function of inverse temper- 
ature /3 in the 3-index matching problem from a population 
dynamics resolution [Tsj of the RS cavity equations @ with 
a large enough value of C. (Here C — 60.). Note that the 
entropy s = I3'^df/df3 is negative for /3 > /3s = 0.412 ± 0.001. 
Inset: overlap q between equilibrium configurations as a 
function of /3; in the glassy phase, the overlap is given by 
q{P,) = 0.321 ±0.002. 



in terms of the probability, exp[/?(xi — C'/d)], of unoccu- 
pation of each of its neighbours (say neighbour i) in the 
Md sites problem. Assuming a replica symmetric (RS) 
structure, the order parameter is the probability V{x) 
that a randomly chosen site i has Xi = x, which satisfies 
the self-consistent equation: 



k d-1 



(5) 



The free energy /rs(/3) can be obtained from 'P{x) as: 



-(In (l + e-'^(«-S.-i-^-))) 



(6) 

where (.) stands for the averages of the cavity fields x 
with the distribution V, of the connectivities k with the 
Poissonian distribution of mean C and of the truncated 
costs ^ with the uniform distribution in [0, C], as in ||SJ). 
In the zero temperature limit, P —^ oo, we obtain formu- 
lae for the ground state energy that directly generalize 
the ones of the BMP case [7llla|. 

However, while correctly describing the d = 2 problem, 
these RS equations are inconsistent when d > 2. We shall 
discuss specifically the d — 3 case. First, the entropy be- 
comes negative for (3 > Ps — 0.412 ± 0.001 , as shown 
on Fig. 13 Second, we have found the RS solution to be 



unstable for /3 > /?,, with Pi ~ 0.6 These two facts 
show that a discontinuous phase transition takes place at 
some inverse temperature Pc ^ Ps- Such transitions are 
also present in other NP-hard combinatorial optimiza- 
tion problems like iiT-SAT, and are usually overcome by 
passing to a one-step replica symmetry broken (IRSB) 
formalism Tl. Here however, the direct application of 
the IRSB cavity method at zero temperatm-e i turns 
out to be inadapted. 

The originality of the MIMP comes from the pecu- 
liar nature of the low temperature phase. This phase 
is dominated by isolated configurations, instead of clus- 
ters of configurations that generally arise in IRSB sys- 
tems Q: the IRSB clusters have no internal entropy 
here, a situation which is also found in some other dis- 
ordered systems, the REM (random ene rgy model [l7||). 
the directed polymer on disordered tree |l8l | and the bi- 
nary perceptron 19]. Upon cooling, these systems freeze 
when reaching the temperature 1/Ps where the entropy 
becomes zero. As a result, the thermodynamical proper- 
ties can be derived from the knowledge of the RS solution 
only |l6|| . The free energy is given by: 



/rs(/3) if /? < 

MPs) itP>Ps. 



(7) 



Necessary conditions for this frozen IRSB Ansatz to 
hold include the existence of a finite Ps where the RS 
entropy becomes negative, the stability of the RS solution 
up to (at least) Ps and the absence of any discontinuous 
IRSB transition before Ps (as we have checked from a 
finite P IRSB population dynamics investigation). 

On top of these properties, a crucial necessary condi- 
tion for the frozen IRSB Ansatz to hold is that the system 
must be subject to a restricted class of constraints, that 
we call hard constraints [l^ . For matching problems, 
hard constraints reflect the requirement to realize perfect 
matchings and basically mean that the occupancy of a hy- 
peredge is uniquely determined by that of its neighbors; 
this is to be contrasted with the situation in coloring 
for instance, where the color of a site is not necessarily 
uniquely prescribed by the colors of its neighbors. No- 
tice that the d = 2 case satisfies all these requirements, 
except for the fact that Ps = oo. 

The prediction IJJ) yields a ground state energy den- 
sity E^/M = frsiPs) = 3.126 ± 0.002 (see Fig. EJ; 
our B&B numerical estimate is compatible with this 
value considering the systematic effects arising from the 
small M used there. When d = 4, we find similarly a 
ground state energy density E^^/M = 7.703 ±0.002 (with 
Ps = 0.135 ± 0.002); here again the B&B estimate we 
obtained is close to this value. 

Overlaps — The cavity method gives access to the 
typical overlap q between equilibrium configurations, de- 
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FIG. 3: Distribution of overlaps between the ground state and 
the first excited state for d = 3 MIMP for M = 8 to M = 22 
(from bottom to top a,t q = .32). 

fined as 

g = ^ E KZZF (8) 

il,---,id 

with (.) and the overline denoting respectively the ther- 
mal and the disorder averages. This overlap can be ex- 
pressed with the cavity method in terms of the order pa- 
rameter Vix). For d = 3, we find q{f3s) = 0.321 ± 0.002. 
Because of the special nature of the frozen IRSB phase 
at /3 > /?s, we expect that, if we take at random two 
configurations among the r lowest energy configurations, 
their overlap will be equal to g(/3s) with probability one 
(for any fixed r, in the large M limit). 

In order to test this prediction, we have generalized the 
B&B method to get numerically the overlap between the 
ground state and the first excited state. Fig.Oshows the 
disorder averaged distribution of the overlap. The data is 
consistent with a distribution becoming peaked at large 
M at an overlap around 0.32, as theoretically predicted 
from the cavity approach. Note also that the overlaps at 
higher values seem to decay to zero: this is exactly the 
prediction of the absence of configurational entropy, i.e., 
a consequence of the freezing scenario which we argued 
arises in this system. 

Another numerical check of the validity of the scenario 
comes from the measurement of the density Af{E) of con- 
figurations as a function of energy. When {E — Eq)/M 
is small, we find that lnJ\f{E) ~ ItiJJ\E) -f{E - Eq), 
with 7(d = 3) ~ 0.405, and 7(d = 4) ~ 0.14. These 
values of 7 agree with the inverse freezing temperature 
(is found in the cavity method. 

Discussion — We have investigated the thermody- 
namics of the d-index matching problem. For d > 3 it 
diff'ers from the (2-index) matching in being NP-hard and 
having a low temperature glassy phase. Physically, in 



the latter case it is much more difficult to find a second 
low energy configuration in the neighborhood of a first 
one. It would be interesting to study this effect further 
along the lines of [i^, |^ . The glassy phase is of a special 
type, distinct from the one found in other recently solved 
NP-complete decision problems, because it has vanishing 
internal entropy. In this respect, the MIMP behaves as 
a REM il2i], freezing into a few configurations. 

We have derived the full phase diagram; we conjecture 
these results to be exact, and the numerical checks which 
we have performed on relatively small systems, through 
an efficient B&B algorithm, are consistent with the pre- 
dictions. It will be extremely interesting to generalize to 
this problem the rigorous mathematical methods devel- 
oped for the BMP. 
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